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1.  Introduction 

The  so-called  “third  generation”  (3G)  model  for  wind-generated 
surface  ocean  gravity  waves  (Komen  et  al.,  1984,1994;  WAMD1 
Group,  1988;  Tolman,  1991;  Booij  et  al.,  1999;  WISE  Group, 
2007)  is  commonly  used  in  engineering  and  operational  forecast¬ 
ing  applications  (e.g.  Jensen  et  al.,  2002;  Bidlot  et  al.,  2002).  In  such 
applications,  the  most  common  and  lowest  order  quantity  pre¬ 
dicted  is  waveheight,  which  is  directly  related  to  the  total  energy, 
m0,  in  a  wave  spectrum,  £(/,  0),  where  £  is  the  spectral  density  of  sea 
surface  elevation  variance  (energy), /is  the  wave  frequency,  and  0 
the  wave  direction:  m0  =  JE(f,0)dfd9.  Increasingly,  attention  is 
being  given  to  higher  order  features  of  the  wave  spectrum.  This 
is  motivated  by  steady  progress  in  accuracy  of  predictions  of  sim¬ 
ple  waveheight,  by  the  resulting  desire  to  examine  in  greater  detail 
the  accuracy  of  E(f,0)  predicted  by  the  models,  and— perhaps  most 
importantly— by  the  inherent  utility  of  higher  order  quantities  for 
specific  applications.  The  present  paper  is  specifically  concerned 
with  the  accuracy  of  predictions  of  the  width,  in  frequency  space, 
of  the  non-directional  spectrum  E(f)  =  /  £(/,  0)d0.  Accurate  predic¬ 
tions  of  the  spectral  narrowness  are  not  only  of  scientific  interest 
but  also  for  engineering  practice,  where  this  property  plays  an 
important  role  in  design  formulas  (Goda,  1985;  Saulnier  et  al.. 
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2011).  Particular  effort  is  made  to  examine  the  impact  of  a  partic¬ 
ular  approximation— the  Discrete  Interaction  Approximation— on 
the  accuracy  of  these  predictions  of  frequency  width.  This  approx¬ 
imation  is  made  in  the  physics  calculations  of  nearly  all  routine  3G 
wave  model  applications  today,  and  is  introduced  below,  subse¬ 
quent  to  the  introduction  of  the  model  physics  in  general. 

The  3G  wave  models  utilize  a  phase-averaged  (spectral) 
description  of  wave  conditions;  the  dependent  variable  being 
solved  for  in  these  models  is  most  often  the  wave  action  density 
/V  =  £/cr(e.g.  Bretherton  and  Garrett,  1968).  Here,  o  is  the  intrinsic 
angular  frequency  where  “angular”  refers  to  the  relation  a  =  2nf. 
Wave  action  N  is  solved  in  five  dimensions,  e.g.  N  =  N(f,0\  x,t), 
where  x  denotes  position  and  t  time.  The  evolution  of  the  wave 
spectrum  is  described  by  means  of  the  radiative  transfer  equation 
(Gelci  et  al.,  1956;  Hasselmann,  1960),  which  can  be  written  as: 

—  \  V  CN  —  —  ^lf1  4  50fher 

dr  ~  o  ~~  a  ^  ' 

where  c  is  the  energy  propagation  velocity  of  the  waves  in  each 
dimension  (c*,cy,ca,Q,).  The  LHS  of  this  equation  accounts  for  kine¬ 
matics,  which  are  conservative,  whereas  on  the  RHS,  Sto[  represents 
dynamics.  In  deep  water,  it  is  generally  accepted  that  wind-wave 
growth  is  primarily  a  result  of  three  physical  processes:  atmo¬ 
spheric  input  from  the  wind  to  the  waves,  Sln,  wave  dissipation 
due  to  breaking  S</5,  and  nonlinear  energy  transfer  between  the 
wave  components  SnI4.  In  finite  depths,  additional  terms  due  to 
wave-bottom  friction,  depth-limited  breaking  and  triad  interactions 
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become  significant,  and  more  terms  can  be  formulated  which  may 
become  important  in  particular  circumstances,  e.g.  non-breaking 
dissipation  or  scattering  from  interaction  with  irregular  bathymetry 
or  with  sea  ice.  Every  3G  model  includes  the  first  three  source 
terms,  whereas  the  remaining  source  terms  can  vary  from  one  mod¬ 
el  to  another.  All  of  these  source  terms  are  spectral  functions.  The 
reader  is  referred  to  WISE  Group  (2007)  for  a  more  complete  over¬ 
view  of  this  technique  of  wave  modeling. 

Of  these  source  terms,  the  present  paper  is  primarily  concerned 
with  the  term  for  four-wave  nonlinear  interactions,  S„f4  (e.g. 
Hasselmann,  1962).  In  a  developing  sea,  this  term  moves  energy 
from  frequencies  just  above  the  spectral  peak  to  both  lower  and 
higher  frequencies.  The  transfer  to  lower  frequencies  is  a  primary 
mechanism  for  frequency  downshifting.  The  transfer  to  higher  fre¬ 
quencies  is  balanced  by  wind  input  and  white-capping  dissipation 
resulting  in  power  law-decay  of  wave  action  with  frequency.  This 
term  has  also  been  found  to  have  a  primary  role  in  determining 
and  stabilizing  the  spectral  shape  in  general  (Hasselmann  et  al. 
1 973 ;  Young  and  Van  Vledder,  1 993).  Unlike  the  other  two  primary 
source  terms,  the  form  of  this  term  is  known  and  can  be  solved  for 
exactly  (e.g.  EXACT-NL,  Hasselmann  and  Hasselmann,  1985)  but  is 
far  too  expensive  for  routine  use,  necessitating  shortcuts,  the  most 
commonly  used  method  being  the  Discrete  Interaction  Approxima¬ 
tion  (D1A)  introduced  in  the  seminal  work  by  Hasselmann  et  al. 
(1985),  which  uses  a  relatively  small  subset  of  interactions  but  pre¬ 
serves  the  lowest  order  features  of  the  exact  approach.  Since  this 
period  of  rapid  progress  in  the  1980s,  much  effort  has  been  made 
in  the  wave  model  development  community  to  alternately  im¬ 
prove  the  accuracy  of  the  DIA  or  improve  the  efficiency  of  the  exact 
algorithms  (WISE  Group,  2007). 

The  qualitative  impact  of  the  DIA  ’’shortcut”  on  model  results 
has  most  often  been  presented  for  idealized  scenarios  (e.g.  fetch- 
limited  or  turning  winds),  as  this  obviously  facilitates  comprehen¬ 
sion  of  the  most  fundamental  features.  Hasselmann  et  al.  (1985), 
for  example,  evaluate  fetch-limited  growth  curves  using  the  DIA 
and  note  that  agreement  of  the  more  important  scale  parameters, 
total  energy  and  peak  frequency,  were  excellent,  while  two  less 
important  scale  parameters  were  predicted  less  well;  the  one  relat¬ 
ing  to  spectral  peakedness  was  reported  as  ’’somewhat  smaller" 
than  the  ground  truth.  In  subsequent  years,  considerable  success 
has  been  achieved  with  3G-modeIs,  especially  in  the  context  of 
those  ’’more  important  scale  parameters",  while  at  the  same  time 
the  shortcomings  of  the  DIA  are  receiving  greater  scrutiny  (e.g. 
Van  Vledder  et  al.,  2000).  These  shortcomings  are  compensated 
by  tuning  of  other  source  terms,  thus  hampering  the  development 
of  3G-models.  In  the  simplest  examples,  the  limitations  of  the  DIA 
are  illustrated  by  computing  Sn,4  using  DIA  for  a  fixed  wave 
spectrum  and  comparing  this  with  exact  computations  of  Sn<4.  This 
exercise  reveals  basic  characteristics  but  is  insufficient  to  antici¬ 
pate  practical  outcomes  in  a  time-stepping  solution,  because  of 
the  highly  non-linear  character  of  SnW.  Instead,  full  model  exercises 
are  needed  in  which  Sn/4  is  used  in  conjunction  with  other  source 
terms.  Some  recent  effort  has  been  made  to  characterize  the 
practical  implications  of  the  shortcut  on  realistic  applications: 
Alves  et  al.  (2002)  examine  fetch-limited  spectra  generated  with 
implementations  of  DIA  and  exact-Snl  in  a  3G  model,  concluding 
that  the  use  of  the  DIA  algorithm  constrains  significantly  the 
shape  of  1  -D  frequency  spectrum;  Tolman  (201 1 )  presents  applica¬ 
tions  in  a  synthetic  hurricane  and  a  storm  event  in  Lake  Michigan. 
One  feature  reported  often  is  broader  directional  spreading  with 
the  DIA  than  with  the  full  solution  (Young  et  al.,  1987;  Van 
Vledder,  1990;  Forristall  and  Ewans,  1998;  Ardhuin  et  al.,  2007; 
Tolman,  2011),  broader  frequency  spectra  (already  visible  in 
Hasselmann  et  al.,  1985,  their  Figure  9)  and  an  underestimation 
of  the  energy  level  at  the  spectral  peak  has  also  been  reported 
(Tolman,  2011). 


Tolman  (201 1 )  and  others  use  the  results  from  ’exact’  calcula¬ 
tions  of  SnM  as  ground  truth.  The  soundness  of  this  approach  is 
obvious,  since  the  accuracy  of  the  ground  truth  cannot  be  ques¬ 
tioned.  However,  these  comparisons  do  not  answer  a  particular, 
interesting  question,  which  is,  "can  evidence  of  these  known  short¬ 
comings  of  the  DIA  be  noted,  isolated,  and  proven  in  comparisons 
to  observational  data?”.  This  is  the  fundamental  question  of  the 
present  paper. 

Rogers  and  Wang  (2007),  henceforth  denoted  ”RW07”  use  a 
multi-month  Lake  Michigan  hindcast  with  the  SWAN  model  (Booij 
et  al.,  1999)  to  investigate  possible  bias  in  predictions  of  directional 
spreading  relative  to  buoy  observations.  Their  study  was  primarily 
motivated  by  previous  claims  that  use  of  the  DIA  leads  to  overpre¬ 
diction  of  directional  spreading.  RW07  demonstrate  overprediction 
of  directional  spreading  above  the  peak  in  idealized  simulations 
(using  the  ’exact’  computations  as  ground  truth)  (their  Figures  2 
and  3),  but  did  not  find  any  overprediction  of  directional  spreading 
in  comparison  to  the  buoy  measurements  (their  Figures  8  and  9).  A 
more  positive  but  less  emphasized  feature  of  RW07  was  the  consis¬ 
tency  of  the  results  with  regard  to/requency  width  (the  topic  of  the 
present  paper).  Specifically,  in  both  the  idealized  comparisons 
against  the  ‘exact’  model  (their  Figure  2)  and  the  buoy  comparisons 
(their  Figures  6  and  7),  there  is  an  overprediction  of  energy  below 
the  spectral  peak,  and  so  the  spectrum  is  too  broad. 

RW07  could  only  indirectly  attribute  overprediction  of  fre¬ 
quency  spreading  to  the  DIA,  since  they  did  not  apply  the  ‘exact’ 
Sn/4  computations  in  their  hindcast  due  to  computational  costs. 
This  is,  however,  done  by  Ardhuin  et  al.  (2007)  in  a  hindcast  for 
Duck,  North  Carolina.  This  was  made  possible  (from  a  practical 
point  of  view)  through  the  use  of  a  finite  element  grid  with  rela¬ 
tively  few  grid  points  (CREST  model,  Ardhuin  et  ah,  2001).  These 
authors  find  that,  in  comparison  to  the  ‘exact’  SnM-based  model, 
the  DIA-based  model  always  has  higher  directional  spreading  (as 
expected),  but  this  does  not  always  lead  to  a  larger  error  in  com¬ 
parison  to  buoy  data.  In  cases  where  directional  spreading  predic¬ 
tions  are  made  worse  with  ’exact’  Sn/4,  this  is  reasonably  attributed 
to  inaccuracies  in  the  other  source  terms.  Ardhuin  et  al.  (2007)  also 
mention  in  passing  that  the  DIA  would  result  in  an  overestimation 
of  the  growth  of  low-frequency  waves  if  not  for  cancellation  of 
errors  via  the  other  source  terms,  thus  leading  to  overly  broad 
frequency  spectra. 

Ardhuin  et  al.  (2010)  conducted  Lake  Michigan  simulations  sim¬ 
ilar  to  RW07,  using  the  WAVEWATCH  III*  model  (”WW3”,  Tolman 
et  al.,  2002;  Tolman,  2009).  These  simulations  indicate  overpredic- 
tion  of  energy  below  the  spectral  peak  similar  to  those  shown  in 
RW07  Figs.  6  and  7.  Again,  this  was  not  positively  attributed  to 
the  DIA.  However,  in  the  context  of  the  Lake  Michigan  hindcasts, 
it  suggests  a  robust  feature,  independent  of  modeling  platform 
(SWAN,  WW3)  and  Sin  +  Sds  parameterization  (Komen  et  al., 
1984;  Bidlot  et  al.,  2005;  Ardhuin  et  al.,  2010,  or,  as  we  will  show, 
Rogers  et  al.,  2012).  The  objective  of  the  present  study  is  to  deter¬ 
mine  whether  this  overprediction  of  frequency  width  is  positively 
attributable  to  the  DIA. 

Direct,  side-by-side  comparisons  of  one-dimensional  wave 
spectra  can  be  a  useful  method  of  presenting  differences,  or  a 
means  of  careful  model  evaluation.  However,  in  the  context  of  rou¬ 
tine  or  repetitive  model  validation  exercises,  it  is  often  impractical. 
This  is  especially  true  in  longer  hindcasts  with  multiple  locations  of 
spectral  observational  data.  In  these  cases,  it  is  useful  to  utilize 
bulk  parameters  that  effectively  identify  model  error,  and  statistics 
can  be  calculated  from  these.  For  example,  zero-moment  wave 
height,  Hm0  is  used  to  quantify  total  energy,  and  the  spectral  mean 
period  (e.g.  Trn._i#0.  Tm>0^.  see  Appendix)  is  used  to  quantify  a 

center  (or  centroid)  of  the  frequency  spectrum,  variously  defined  to 
give  more  or  less  weight  to  frequencies  further  from  that  center.  In 
this  paper,  existing  methods  for  quantifying  the  narrowness  in 
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frequency  space  of  non-directional  frequency  spectra  (as  opposed  to 
spectral  width)  are  evaluated.  We  find  that  some  methods  are  not 
sufficiently  sensitive  to  problems  with  narrowness  that  are  readily 
visible  by  eye. 

This  paper  is  organized  as  follows.  The  wave  model  is  described 
in  Section  2,  and  the  utilized  hindcast  is  described  in  Section  3.  The 
methods  of  quantifying  frequency  narrowness  are  described  in 
Section  4.  In  Section  5,  overprediction  of  spectral  width  is  demon¬ 
strated  in  a  hindcast  using  DIA.  In  Section  6,  an  ‘exact’  Sn/4-based 
hindcast  is  analyzed  and  it  is  found  that  the  bias  in  frequency 
width  is  reduced.  Discussion  and  Conclusions  are  given  in  Sections 
7  and  8. 

2.  Mode!  description 

The  model  platform  is  SWAN  (Booij  et  al.t  1999;  SWAN,  2010) 
and  the  wind  input  (S,*n),  whitecapping  (5*),  and  non-breaking  dis¬ 
sipation  (Ssweu)  parameterizations  are  taken  from  Rogers  et  al. 
(2012).  We  refer  the  reader  to  that  paper  for  details,  but  the  key 
features  are  briefly  given  here.  The  wind  input  term  S,„  is  based 
on  Donelan  et  al.  (2006)  and  Tsagareli  et  al.  (2010),  developed  from 
direct  measurements  of  wind  input  at  Lake  George,  Australia,  with 
the  wind  drag  coefficient  Cd  based  on  Hwang  (2011).  Dissipation 
from  breaking  (whitecapping),  Sd$  is  based  on  Babanin  et  al. 
(2010)  which  is  developed  from  Young  and  Babanin  (2006), 
Tsagareli  (2009),  and  Banner  et  al.  (2000).  Srfs  is  two-phase— insofar 
as  it  accounts  for  breaking  of  waves  due  to  inherent  instability  and 
dissipation  induced  by  the  breaking  of  longer  waves— and  employs 
a  breaking  threshold,  based  on  the  Phillips  (1958)  saturation  spec¬ 
trum  (the  same  concept  was  used  earlier  by  others,  such  as  Collins 
(1972)  and  Alves  and  Banner  (2003)).  In  Rogers  et  al.  (2012),  the 
utilized  form  of  S^s  is  denoted  as  L4M4.  Non-breaking  dissipation, 
Sswcii .  is  included  to  account  for  the  slow  attenuation  of  swell  by 
non-breaking  processes,  and  utilizes  work  by  Ardhuin  et  al. 
(2009,2010).  In  the  notation  of  Rogers  etal.  (201 2 ),  fe  =  0.006  (con¬ 
trolling  the  strength  of  swell  dissipation)  is  used  in  simulations  de¬ 
scribed  below.  For  four-wave  nonlinear  interactions  Sn(4,  we  use 
methods  available  in  SWAN  (2010)  without  modification.  The  term 
is  computed  with  either  DIA  or  ‘exact’  computations.  The  latter  is 
denoted  “XNL”  herein  and  is  specifically  the  “Webb-Resio-Tracy” 
method  (WRT)  as  implemented  by  van  Vledder  in  SWAN,  see  Van 
Vledder  (2002,  2006)  and  SWAN  (2010). 


3.  Hindcast  description:  Lake  Michigan  2002, 10-day  simulation 

Lake  Michigan  is  a  useful  area  for  studying  the  impact  of  model 
physics,  as  it  is  mostly  deep  water  (minimizing  complications  due 
to  finite  depth  source  terms),  large  enough  to  allow  generation  of 
dominant  waves  measurable  using  National  Data  Buoy  Center 
(NDBC)  buoys  (e.g.  Tp  =  6s),  and  enclosed,  minimizing  complica¬ 
tions  due  to  interaction  with  swell  (older  swells  being  non-exis¬ 
tent).  The  latter  feature  is  pertinent  to  the  scope  of  this  study, 
which  is  specific  to  windsea. 

The  Lake  Michigan  model  simulations  are  set  up  as  follows.  The 
computational  grid  is  approximately  252  km  (east-west)  by 
496  km  (north-south),  with  4  km  grid  resolution.  Directional  reso¬ 
lution  is  10°  (36  bins),  and  a  logarithmic  frequency  grid  is  used, 
with  35  frequencies  from  0.07  to  1.97  Hz.  The  bathymetry,  shown 
in  Fig.  1,  is  provided  at  2  km  resolution  by  the  NOAA  Great  Lakes 
Environmental  Research  Laboratory.  Wind  forcing  is  created  from 
two  NOAA  NDBC  buoys  as  described  in  Rogers  et  al.  (2003)  and 
RW07;  the  method  assumes  homogeneity  in  longitude  and  spatial 
variation  in  latitude  determined  by  the  buoys  45002  and  45007. 
JONSWAP  (Hasselmann  et  al.  1973)  bottom  friction  formula  is 
used,  though  this  is  not  expected  to  affect  results,  since  this 


Fig.  1.  Bathymetry  of  Lake  Michigan  (depths  in  meters).  NOAA  NDBC  station 
locations  are  indicated.  Buoy  45007  is  used  for  model/data  comparisons  herein. 


hindcast  is  predominately  in  deep  water.  The  default  numerics 
are  used,  a  second  order  scheme  with  third  order  diffusion  (see 
SWAN,  2010).  Physics  for  Sin,  Sds,  and  S^u  are  used  as  described 
in  Section  2.  The  initial  condition  is  from  a  ’hotstart’  file  created 
by  the  model  using  DIA  for  SnM,  0000  UTC  12  Oct.  2002  to  0000 
UTC  13  Oct.  2002.  The  simulations  used  for  comparisons  are  from 
0000  UTC  13  Oct.  2002  to  0000  UTC  23  Oct.  2002  (10-day  dura¬ 
tion),  using  a  time  step  of  6  min. 


4.  Methods 


As  mentioned  in  Section  1,  the  objective  of  the  present  study  is 
to  determine  whether  overprediction  of  frequency  width  is  posi¬ 
tively  attributable  to  the  DIA.  A  necessary  step  is  to  evaluate  the 
utility  of  bulk  parameters  which  can  be  used  to  quantify  the  nar¬ 
rowness  in  frequency  space  of  the  non-directional  spectrum.  In 
this  section,  four  methods  are  described  to  quantify  this  property. 
All  four  methods  are  based  on  simple,  direct  calculation,  mostly  via 
integration,  which  is  preferred  over  indirect  methods,  e.g.  fitting  to 
parametric  spectra;  this  is  a  subjective  decision  and  we  do  not  dis¬ 
pute  the  merits  of  methods  of  the  latter  type. 

Method  A:  Uses  the  method  of  SWAN  (2010),  as  defined  by  Batt- 
jes  and  Van  Vledder  (1984).  The  calculation  is: 


ffmox 

Qa  =  /  E(f)ea'x,xdf 

I  ”  fmh i 


/m0 


(2) 


where  m0  =  ff™  E(J)df  and  r  =  Tm,0^  (defined  in  Appendix).  The  low¬ 
er  and  upper  bounds  of  the  integration  are  denoted  as/mm  and  fmax. 
As  with  other  integral  quantities,  for  any  given  comparison,  fmin  and 
fmax  should  be  applied  consistently:  when  comparing  model  and 
observations,  this  will  often  imply  that  the  highest  model  frequen¬ 
cies  are  excluded  from  the  calculation.  In  Battjes  and  Van  Vledder 
(1984),  Q*  is  denoted  as  parameter  k  and  is  referred  to  as  a  “shape 
parameter”  for  the  purpose  of  predicting  wave  group  statistics  from 
a  non-directional  spectrum.  The  Battjes  and  Van  Vledder  (1984) 
paper  is  given  as  the  source  of  the  equation  in  SWAN  (2010). 
although  the  actual  source  is  Rice  (1944).  However,  in  SWAN 
(201 0),  the  quantity  is  incorrectly  referred  to  as  “FSPR”,  “the  normal¬ 
ized  frequency  width  of  the  spectrum  (frequency  spreading)”. 
Like  the  other  methods  described  here,  this  parameter  actually 
quantifies  the  narrowness  of  the  spectrum  and  it  ranges  from  0  to  1. 
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Method  B  is: 

Qb  =  (^m.O,2/Tmt-1.o)  (3,» 

where  Tm0j  and  Tm._it0  are  defined  in  the  Appendix.  Thus,  it  is  the 
square  of  the  ratio  between  Tm,0.2.  which  is  often  70-90%  of  the  peak 
period,  and  Tm_lt0f  which  is  usually  much  closer  to  the  peak  period 
The  closer  this  parameter  is  to  one,  the  more  peaked,  or  narrow,  the 
frequency  spectrum. 

A  similar  quantity,  Tmto.i/Tmt.  it0  was  used  by  Van  Vledder  et  al. 
(2008),  in  studies  of  the  Dutch  Waddensea  to  assist  in  detection 
of  areas  where  longer  waves  were  dominant  and  where  bi-model 
spectra  occurred.  In  the  numerator  of  Qb,  Tm,0.2  is  used  here  in  pref¬ 
erence  to  Tm> o,i  as  it  provides  more  separation  from  and 

thus  it  is  expected  to  have  greater  sensitivity  to  changes  in  shape. 
Method  C  is: 


Qc  =  max(E/„)/Tm  _,i0  (4) 

where: 


Efn(f)  =  E(f)/m0  (5) 

and: 


rfmax 

/  E/n(f)df  = 

jfmin 


(6) 


Thus  Qc  is  the  peak  of  the  function  E(J)  after  it  has  been  normalized 
by  the  area  under  the  same  function.1  The  division  by  Tm>  lj0  is 
included  to  make  the  quantity  dimensionless.  This  method  was 
devised  in  this  study  following  the  formulation  of  the  A  parameter 
used  by  Babanin  and  Soloviev  (1987,  1998a)  to  quantify  narrowness 
of  spectra  in  directional  space.  However,  with  the  non-dimensional- 
ization,  we  find  that  the  formula  is  similar  to  those  used  by  LeBlond 
et  al.  (1982),  Belberov  et  al.  (1983),  Mansard  and  Funke  (1990),  and 
Babanin  and  Soloviev  (1998b).  Of  the  four  methods  used  here,  only 
this  method  includes  a  non-integral  operation,  the  Mmax(E/n)"  opera¬ 
tion  in  (4)  This  implies  that  Qc  responds  to  spectral  peakedness. 

Method  D  is: 


(b^2/m20  [U'’fE2(J)df  (7) 

V/mm 

This  is  taken  from  Goda  (1970,  1985).  It  has  been  referenced  in  the 
freak  wave  literature  (Janssen,  2003)  as  "Coda’s  peakedness  factor 
Qp”  with  the  Benjamin  Feir  Index  being  BFI  =  k0Qp y/rn^In.  The  high¬ 
er  Qp,  the  more  peaked  the  spectrum. 

As  it  is  less  useful  to  evaluate  higher  order  moments  (in  this 
case,  spreading)  in  cases  where  the  lower  order  moments  are  dis¬ 
similar,  we  only  include  points  in  the  time  series  for  which  eH  <  0.2 
and  er<  0.2,  where: 

“  ^mO.obs  ~  ^mO, model] /^mO.obs 
Cr  =  (T'm.-  1 .0.o  bs  —  7'm,-1.0.modt/|/7'm.-1.0.obs 

To  familiarize  with  these  quantities,  we  compare  in  Table  1  values 
calculated  using  the  JONSWAP  spectrum  (see  for  example.  Young, 
1999,  pg.  112)  for  different  values  of  the  JONSWAP  spectral  peaked¬ 
ness  parameter  y,  which  many  readers  will  be  familiar  with.2  JON¬ 
SWAP  y*l  corresponds  to  fully  developed  seas,  which  tend  to  be 
relatively  broad  in  frequency  space,  and  JONSWAP  7 -3.3  corre- 


1  Apart  from  the  division  by  the  spectral  wave  period,  one  can  recover  the  same 
value  using  the  reciprocal  value  of  the  area  of  the  function  that  is  normalized  by  the 
peak  value. 

2  Though  it  is  not  expected  to  have  much  impact  on  the  calculations,  for 

completeness,  the  other  JONSWAP  parameters  used  are:  /p-0.19  and  3-0.02. 
cram  0.07,  -  0.09  (see  e.g..  Young,  1999  for  definition)  with  integration  of 

frequencies  from  fmSn  -0.03  Hz  with  grid  spacing  A /•  0.001  Hz  and  an  /  5  spectral 
tail. 


Table  1 

Example  calculations  of  frequency  narrowness  parameters  on  JONSWAP  spectra.  In 
each  case,  two  values  for  Q  are  given.  The  first  value  corresponds  to  upper  bound  of 
integration  fmaMm  1.0  Hz  and  the  second  value,  in  round  brackets  (),  corresponds  to 
fmax  -  0.4  Hz. 


JONSWAP  j 

?  Narrowness 

Qa 

Qb 

Qc 

Qd 

1.0 

0.42  (0.37) 

0.72  (0.83) 

1.67  (1.72) 

2.01  (2.25) 

2.0 

0.50  (0.47) 

0.74  (0.85) 

2.60  (2.66) 

2  49  (2.74) 

3.3 

0.57  (0.55) 

0  77  (0.87) 

3.43  (3.49) 

3.15(3.41) 

6.0 

0.65  (0.65) 

0.80  (0.89) 

4.54  (4.60) 

4.32  (4.59) 

10.0 

0.72  (0.72) 

0.84(0.91) 

S.54  (5.S9) 

5.62  (5.87) 

100.0 

0.92  (0.92) 

0.96  (0.98) 

9.56  (9.58) 

12.15  (12.25) 
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Fig.  2.  Comparison  of  observed  and  computed  non -directional  spectra  using  the 
simulation  with  DIA  for  Sn/4  for  six  time  periods. 


sponds  to  younger  seas,  which  tend  to  be  relatively  narrow  in  fre¬ 
quency  space.  Values  of  y  >  6  are  not  necessarily  realistic  for  wind 
seas,  but  are  included  here  for  illustration.3  Also  note  that  peak 
enhancement  is  a  specific  quantity  associated  with  "overshoot"  near 
the  spectral  peak  (e.g.  Young  (1999),  Figure  5.20),  whereas  spectral 
narrowness  is  a  more  general  quantity  that  might  be  affected  by 
other  features  of  the  spectrum,  e.g.  a  bimodal  frequency  spectrum 
will  tend  to  have  a  low  spectral  narrowness. 

All  four  quantities  are  dimensionless.  Separate  calculations  with 
the  JONSWAP  spectra  (not  shown)  indicate  that  none  of  the  Q 
parameters  exhibit  a  significant  dependence  on  the  peak  frequency 
fp.  With  frequency  distributions  shaped  like  normal  distributions, 
Qc  and  Qp  exhibit  strong  and  identical  dependence  on  fp  (linear 
proportionality),  a  logical  consequence  of  the  form  of  the  equa¬ 
tions.  Qa  and  Qb  also  exhibit  some  dependence  on  fp  using  normal 
distributions.  The  differing  behavior  associated  with  use  of  normal 
distributions  versus  JONSWAP  spectra  is  caused  by  the  dependence 
of  the  JONSWAP  spectrum  on/p  beyond  simply  shifting  the  spec¬ 
trum  in  frequency  space. 

The  narrowness  values  in  Table  1  are  calculated  with  two  alter¬ 
nate  values  for  the  upper  limit  of  integration,  /max«  1.0  Hz  and 
fmax  m  0.4  Hz.  This  reveals  some  dependency  of  the  parameters  on 
fmax ♦  Method  B  in  particular  is  sensitive,  which  is  not  surprising  gi¬ 
ven  that  it  utilizes  a  higher  order  moment  of  the  spectrum.4  Meth- 


3  Using  y  >  6  is  actually  a  technique  for  providing  SWAN  with  parameterized 
boundary  forcing  that  is  narrow,  i.e.  an  unrealistically  peaked  wind  sea  can  be  used  to 
approximate  older  swells,  since  the  latter  tend  to  be  narrow  in  the  real  ocean. 

4  The  sensitivity  would  presumably  be  reduced  by  using  Tm>0.i  instead  of  Tm0J. 
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Fig.  3.  Energy  prediction  estimates  (4v/m3,  in  meters),  divided  into  four  relative 
frequency  bands:  DIA-based  model  vs.  observations.  Statistical  quantities  shown 
are  in  order:  number  of  observations,  bias,  root-mean-square  error,  normalized 
root-mean-square  error  (calculation  given  in  the  Appendix),  scatter  index  (standard 
deviation  of  errors  divided  by  the  mean  of  observations),  and  correlation  coefficient 
CC  calculated  as  shown  in  the  Appendix. 


buoy 

frequency  narrowness  Q 


buoy 

frequency  narrowness  Q 


Fig.  4.  Validation  of  frequency  narrowness  parameters  Q*.  Qb.  Or.  and  Qo-  All  model 
results  shown  are  from  the  DIA  simulation.  The  red  line  shows  the  linear  best-fit 
relation.  Statistical  quantities  shown  are  described  in  the  Appendix.  (For  interpre¬ 
tation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the 
web  version  of  this  article.) 


od  C  appears  to  be  the  least  sensitive.  For  example,  the  relative 
change  from  fmax  =  1.0  Hz  to  fmax  =  0.4  Hz  with  JONSWAP  y  =1.0  is 
12%,  +15%,  +3%,  and  +12%,  for  Q^,  Qb,  Qc  and  Q^,  respectively.  At 
y  -  6.0,  it  is  0%,  +11%,  +1%,  and  +6%. 

5.  Results  using  DIA 

In  this  section,  results  from  the  model  which  employs  the  DIA 
for  Sn/4  are  compared  to  observations  from  NOAA  NDBC  buoy 
45007.  Traditional  bulk  parameters  such  as  waveheight  and  wave 
periods  are  only  of  peripheral  concern  in  this  study,  but  for  sake 
of  completeness,  some  are  evaluated  in  the  Appendix. 

As  a  first  step,  the  most  direct  method  of  comparison  is  used. 
Non-directional  spectra  E(f)  are  compared  for  six  time  periods  dur¬ 
ing  the  10-day  simulations  in  Fig.  2.  These  six  times  are  selected  as 
periods  in  which  the  absolute  error  of  Q^, 


^Qc  “  Qc.model  Qc, 


\0b5f 


(9) 


matches  the  median  Eqc  for  the  10-day  simulation.  Thus,  this  small 
subset  of  the  larger  simulation  can  be  regarded  as  a  typical  repre¬ 
sentation  of  the  whole,  in  the  context  of  the  narrowness  quantity 
Qc.  Features  of  the  DIA-based  simulation  are  noticeably  similar  to 
those  noted  by  RW07,  Ardhuin  et  al.  (2007)  and  Tolman  (201 1 ).  Rel¬ 
ative  to  the  observations:  (1)  the  energy  level  at  the  peak  is  low  in 
five  of  the  six  plots,  and  (2)  there  is  too  much  energy  below  the  peak 
in  at  least  three  plots.  However,  as  mentioned  in  Section  1,  such 
qualitative  comparisons  cannot  be  directly  used  to  generate  statis¬ 
tics  for  long  time  series.  Another  form  of  comparison  is  to  look  at 
energy  level,  separated  according  to  frequency  relative  to  the  peak, 
such  as  done  in  RW07  (their  Figure  6).  We  can  thereby  verify  that 
the  DIA-based  model  exhibits  the  same  overprediction  of  energy 
below  the  peak,  as  observed  by  that  earlier  study.  This  is  shown 
in  Fig.  3;  the  quantity  compared  is: 


HmOjxmial  — 


(10) 


where  frequencies/!  and/2  are  indicated  in  each  panel.  The  bias  in 
the  frequency  range  of  0.5  fp  </<  0.8  fp  is  +14  cm.5 

The  frequency  narrowness  parameters  introduced  in  Section  4 
are  applied  to  evaluate  the  sensitivity  of  these  four  parameters  to 
the  problems  visible  in  comparisons  such  as  Fig.  2.  This  is  done 
for  the  DIA-based  simulation  results  in  Fig.  4.  In  all  cases,  the  upper 
bound  on  integration,  fmaxi  is  0.4  Hz,  corresponding  to  the  buoy’s 
upper  limit.  We  find  that  parameters  and  Qb  do  not  show  (or 
only  weakly  show)  a  bias  that  is  apparent  in  the  non-directional 
spectral  comparisons.  Parameters  Qc  and  Qd  do  demonstrate  this 
sensitivity.  Qc  indicates  the  largest  normalized  bias,  and  it  was 
demonstrated  in  Section  4  that  this  parameter  has  less  sensitivity 
to  the  selection  of fmax,  so  it  is  selected  for  all  subsequent  compar¬ 
isons,  though  it  is  noted  that  Qo  is  also  suitable. 


6.  Results  using  XNL 

The  10-day  time  simulation  is  repeated  using  XNL  for  SntA.  For 
the  other  source  terms,  no  changes  are  made,  e.g.  there  is  no  re-cal¬ 
ibration  of  the  coefficients  used  in  S*.  Fig.  5  shows  results  from  6  h 
during  the  10-day  simulations,  for  comparison  with  Fig.  2.  Qualita¬ 
tively,  the  XNL-based  simulation  appears  to  provide  a  better  match 
to  the  observations— especially  in  the  context  of  the  energy  level  at 
the  peak,  and  amount  of  energy  below  the  peak— in  at  least  four  of 
the  six  examples.  A  comparison  of  Hm 0pflrtia,,  similar  to  Fig.  3,  is 
made  in  Fig.  6.  Whereas  the  bias  in  the  frequency  range  of 
0-5/p  </<  0.8/p  is  +14cm  using  the  DIA-based  model,  corre¬ 
sponding  results  with  the  XNL-based  model  indicate  much  smaller 
overall  bias  in  energy  below  the  peak, +5  cm.  Thus,  our  data  clearly 
support  the  interpretation  that  inaccuracies  with  the  DIA  contrib¬ 
ute  to  the  overprediction  of  energy  below  the  peak. 

To  provide  further  support  to  this  interpretation,  the  frequency 
narrowness  parameters  introduced  in  Section  4  are  applied.  Figs.  7 
and  8  compare  results  using  XNL  versus  those  using  DIA  in  terms  of 
frequency  narrowness  parameters  Qc  and  Qd.  As  shown  in  Fig.  7, 

5  Using  a  similar  simulation  with  the  Komen  et  a1„  1984  physics,  we  find  a  bias  of 
+1 1  cm;  using  the  74-day  simulation  with  Komen  et  at,  1984  physics,  RW07  report  a 
bias  of  +9  cm. 
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Fig.  5.  Comparison  of  observed  and  computed  non-directional  spectra  using  the 
simulations  with  D1A  and  XNL  for  SnM  for  six  time  periods. 
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Fig.  6.  Similar  to  Fig.  3,  except  for  XNl-based  model. 


the  bias  in  Qc  (upper  panels  of  Fig.  7 )  is  almost  completely  removed 
using  XNL  and  the  best-fit  slope  is  near  unity.  However,  random  er¬ 
ror,  as  quantified  by  the  scatter  index  and  correlation  coefficient, 
actually  increases;  the  cause  of  this  is  not  known.  Both  the  positive 
and  negative  consequences  of  replacing  DIA  with  XNL  are  seen  also 
in  the  Qo  comparison  (lower  panels  of  Fig.  7),  but  the  consequences 
are  not  as  pronounced  as  with  Qc.  The  time  series  comparison  of 
Qc,  shown  in  Fig.  8,  clearly  illustrates  the  improved  results  by  using 
XNL  instead  of  the  DIA, 

7.  Discussion 

Model  performance  is  often  expressed  in  terms  of  bulk  statistics 
of  significant  wave  height  Hm0  and  peak  period  Tp  or  a  mean  spec¬ 
tral  period,  e.g.  Tm#0.2.  For  many  applications  these  first  order  met¬ 
rics  are  sufficient  to  judge  model  performance.  However,  with 
increasing  requirements  on  model  skill,  other  integral  spectral 
parameters  are  gradually  included  in  more  model  verification 


Fig.  7.  Validation  of  frequency  narrowness  using  the  DIA  and  XNL  for  evaluating  the 
quadruplet  interactions.  The  frequency  narrowness  is  quantified  using  parameter 
Qc  (upper  panels)  and  Qo  (lower  panels).  Statistical  quantities  shown  are  described 
in  the  Appendix. 


model  =  v55UL4M4H  fe  0.0060 


MM/DD 


Fig.  8.  As  in  upper  panels  of  Fig.  7,  except  as  time  series  comparison. 


studies,  such  as  the  spectral  period  the  mean  wave  direc¬ 

tion  0m  and  the  circular  RMS  directional  spreading  crfl.  Together 
these  parameters  form  a  comprehensive  set  to  assess  model  per¬ 
formance,  but  their  error  behavior  can  also  be  used  to  pinpoint 
model  deficiencies. 

In  this  study  frequency  narrowness  is  added  to  this  list  to  inves¬ 
tigate  the  effect  of  computational  methods  for  nonlinear  quadru¬ 
plet  interactions  on  the  narrowness  of  the  frequency  spectrum. 
To  that  end  various  metrics  quantifying  spectral  narrowness  were 
investigated  and  one  parameter  was  selected  that  is  able  to  show  a 
specific  deficiency  of  the  Discrete  Interaction  Approximation  on 
the  narrowness  (or  peakedness)  of  the  frequency  spectrum,  and 
the  over-prediction  of  spectral  density  below  the  peak  frequency 
in  particular.  Based  on  simulations  in  Lake  Michigan  using  the 
SWAN  model  and  DIA  and  XNL  parameterization  of  the  quadruplet 


58 


W.E.  Rogers,  G.Ph.  Van  Vledder/ Ocean  Modeling  70  (2013)  52-6? 


Table  2 

Statistics  for  simulations  of  duration  10  days.  The  number  of  spectra  used  is  indicated  as  n.  Columns  2  and  4  correspond  to  the  simulations  discussed  in  the  main  text,  and  the 
physics  described  in  Section  2  are  denoted  as  R8W2012. 


XNL  R8W2012 

D1A  RBW2012 

D1A  KHH  nk-  2 

Full  set 

5ubset 

Full  set 

Subset 

Full  set 

5ubset 

n 

240 

156 

240 

166 

240 

172 

Hm<0  (m) 

8ias 

-0,05 

0.01 

-0.03 

-0.00 

0.00 

0.03 

RMSE 

0.21 

0.13 

0.18 

0.12 

0.17 

0.12 

CC 

0.93 

0.97 

0.95 

0.98 

0.96 

098 

Tn.  0.1  (s) 

8ias 

-0.18 

-0.15 

-0.10 

-0.06 

0.08 

0.15 

RMSE 

0.37 

0.28 

0.33 

0.23 

0.32 

0.26 

CC 

0.85 

0.92 

0.86 

0.94 

0.89 

0.95 

Tn.0.2  (S) 

8ias 

-0.18 

-0.15 

-0.1 1 

-0.07 

0.06 

0  12 

RMSE 

0.35 

0.27 

0.31 

0.21 

0.29 

0.23 

CC 

0.85 

0.91 

0.86 

0.94 

0.89 

0.95 

Im.  i,o  (S) 

8ias 

-0.18 

-0.14 

0.07 

-0.01 

0.12 

0.20 

RM5E 

0.41 

0.31 

0.37 

0.26 

0.38 

0.33 

CC 

0.85 

0.92 

0.86 

0.93 

0.88 

0.94 

0* 

Norm,  bias 

-0.04 

-0.06 

0.14 

-0.16 

-0.09 

-0.11 

CC 

0.43 

0.51 

0.42 

0.43 

0.46 

0.44 

Slope 

0.94 

0.94 

0.86 

0.86 

0.91 

0.90 

Qb 

Norm,  bias 

-0.00 

-0.01 

-0.02 

-0.03 

-0.02 

-0.03 

CC 

0.74 

0.79 

0.69 

0.72 

0.69 

0.72 

51ope 

1.00 

0.99 

0.98 

0.97 

0.98 

0.97 

Qc 

Norm,  bias 

-0.01 

0.01 

-0.27 

-0.28 

-0.26 

-0.26 

CC 

0.34 

0.38 

0  47 

0.52 

0.49 

0.49 

Slope 

0.96 

0.98 

0.71 

071 

073 

0.73 

Qd 

Norm,  bias 

-0.01 

0.03 

-0.17 

-020 

-0.16 

-0.18 

CC 

0.53 

0.59 

056 

0.66 

059 

0.65 

5lope 

095 

0.95 

0.80 

0.79 

0.82 

0.81 

interactions,  we  are  able  to  show  the  applicability  of  the  spectral 
narrowness  parameter  Ck  to  illustrate  a  potential  weak  point  of 
the  DIA* 

Other  metrics  for  frequency  width,  not  shown  here,  have  been 
used  with  success.  For  example,  in  cases  of  uni-modal  spectra, 
the  distance  (in  Hz)  between  fo  and  f2  can  be  used,  where 
H(fi)  -  E(f2)  -  O.Smax  (£(/)).  Or,  the  spectrum  can  be  normalized 
and  treated  as  a  pseudo-probability  distribution  function,  and 
standard  deviation-like  parameter  can  then  be  calculated.  Such 
parameters  might  in  fact  be  preferred  to  the  Q  parameters  since 
the  units  (Hz)  are  more  tangible.  Saulnier  et  al.  (2011)  recently  pro¬ 
vided  an  excellent  overview  of  a  number  of  metrics,  both  dimen¬ 
sional  and  non-dimensional.  As  mentioned  already,  fitting  a  y 
parameter  to  the  JONSWAP  has  been  used  successfully  in  the  past, 
e.g.  by  Hasselmann  et  al.  (1985).  Experiments  with  16  alternate 
metrics  have  been  performed  and  will  be  documented  separately 
in  a  report  focused  on  metrics. 

Validation  of  higher  order  moments  of  directional  and  non- 
directional  spectra  introduces  new  challenges.  Not  only  are  the 
models  normally  less  reliable  for  prediction  of  these  moments, 
but  such  applications  also  require  more  fidelity  from  our  observa¬ 
tional  datasets,  perhaps  pushing  to  their  limits  in  some  cases.  Van 
Vledder  and  Battjes  (1992)  raise  concerns  about  the  statistical 
properties  of  narrowness  metric  Qd,  which  is  based  on  Goda 
(1970,  1985).  In  short,  though  the  metric  can  be  quite  useful  for 
application  to  model  spectra,  application  to  measured  spectra  is 
less  reliable  due  to  sensitivity  to  the  amount  of  smoothing  used 
in  creating  the  spectra.  We  anticipate  that  this  criticism  would  also 


apply  to  method  C  as  it  may  affect  the  peak  value  max(Efn).  As  such, 
validation  with  these  metrics  (such  as  our  Fig.  7)  must  be  inter¬ 
preted  in  the  context  of  the  dataset  used.  Further,  disparate  obser¬ 
vational  datasets  should  not  be  combined  when  creating  statistics 
for  validation:  the  statistics  should  be  treated  separately.  This  is  in 
contrast  to  treatment  of  traditional  quantities  such  as  wave  height. 

Though  much  more  useful  than  subjective  visual  comparison  of 
non-directional  spectra,  comparisons  such  as  Fig.  3  still  have  limita¬ 
tions.  The  most  significant  limitation  is  the  requirement  to  identify 
the  spectral  peak,  and  that  all  energy  is  evaluated  according  to  its  po¬ 
sition  relative  to  this  single  peak.  This  is  easily  justified  in  environ¬ 
ments  such  as  Lake  Michigan,  which  is  typically  dominated  by 
windsea.6  However,  in  the  open  ocean,  it  is  much  more  common  to 
have  multiple  peaks  due  to  swell,  and  so  these  comparisons  are  impos¬ 
sible  without  first  separating  the  windsea  from  swell.  The  four  nar¬ 
rowness  parameters  applied  here,  by  contrast,  still  have  meaning  in 
the  case  of  bimodal  spectra,  just  as  waveheightHm0  has  meaning.  Even 
so.  it  is  felt  that  application  of  a  windsea/swell  separation  algorithm 
(e.g.  Gerling,  1992;  Hanson  and  Phillips,  2001)  prior  to  calculations 
of  narrowness  will  often  be  worthwhile,  and  it  obviously  necessary 
if  windsea  (e.g.  impact  of  source  terms  on  windsea)  is  being  studied. 

In  Section  5,  the  overprediction  of  wave  energy  below  the  spec¬ 
tral  peak  using  a  DlA-based  model  is  documented,  using  the  phys- 


6  This  dominance  is  not  absolute,  however:  in  the  Great  Lakes,  because  of  their  size, 
buoy  measurements  sometimes  indicate  two  peaks  during  rapidly  rotating  winds, 
implying  a  new  windsea  and  old  windsea  (alternately  named  young  swell) 
component 
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Table  3 

Statistics  for  simulations  of  duration  74  days  using  D1A.  The  third  column  (KHH  nk  -  2,  full  set)  is  essentially  the  same  simulation  as  used  by  RW07. 


RBW2012 

KHH  nk  -  2 

KHH  n*-  1 

Full  set 

Subset 

Full  set 

Subset 

Full  set 

Subset 

n 

1685 

1092 

1688 

1089 

1688 

912 

Hm.0  (m) 

Bias 

-0.02 

-0.01 

0.02 

0.01 

-0.12 

-  0.06 

RMSE 

0.17 

0.12 

0.17 

0.12 

021 

0.14 

CC 

0.96 

0.97 

0.96 

098 

0.95 

0.98 

TrUU  (s) 

Bias 

-0.07 

-0.07 

0.11 

0.12 

-0.36 

-0.26 

RMSE 

0.33 

0.24 

0.34 

0.26 

0.51 

0.35 

CC 

0,90 

0.94 

0.90 

0.94 

0.88 

0.95 

Tn.0.2  (S) 

Bias 

0.08 

-0.08 

0.09 

0.10 

-0.35 

-026 

RMSE 

0.31 

0.22 

0.32 

0.24 

0.49 

0.34 

CC 

090 

0.94 

0.90 

0.94 

0.88 

0.95 

Tn.  t.0  (S) 

8ias 

0.05 

-0.03 

0.15 

0.16 

-0.37 

-023 

RMSE 

0.36 

0.26 

0.40 

0.32 

0.55 

0.35 

CC 

0.90 

0.94 

0.91 

0.94 

0.87 

0.95 

Qa 

Norm,  bias 

-0.11 

-0.13 

-0.08 

0.08 

0.08 

-0.14 

CC 

0,76 

0.70 

0.76 

0.68 

0.63 

0.61 

5  lope 

0.88 

0.87 

0.90 

0.92 

0.88 

0.83 

Qfl 

Norm,  bias 

-0.02 

-0.02 

0.02 

-0.02 

-0.00 

-0.02 

CC 

0.84 

0.85 

0.83 

0.84 

0.74 

0.80 

51ope 

0.98 

0.98 

0.98 

098 

0.99 

0.98 

Or 

Norm,  bias 

-0.29 

-028 

-0.32 

-0.25 

-0.30 

-0.30 

CC 

0.53 

0.61 

0.57 

0.63 

0.58 

0.55 

51ope 

0.70 

0.70 

0.71 

0.73 

0.69 

0.68 

Qjd 

Norm,  bias 

-0.28 

-0.20 

-0.34 

-0.17 

-027 

-0.19 

CC 

0.46 

0.78 

0.52 

0.81 

0.53 

0.74 

5lope 

0.80 

0.80 

0.80 

0.81 

0.83 

0.78 

ics  of  Rogers  et  al.  (2012)  for  S,-n  and  S*.  The  same  behavior  is  doc¬ 
umented  for  longer  duration  simulations  and  more  traditional 
forms  of  Sin  and  Sds  in  the  Appendix.  Further,  the  behavior  is  ob¬ 
served  by  Ardhuin  et  al,  (2010)  using  a  third  set  of  physics,  and  a 
different  model,  WAVEWATCH  HI,  applied  in  Lake  Michigan.  Ard¬ 
huin  et  al.  (2007)  observe  similar  behavior  using  a  third  model,  ap¬ 
plied  at  the  North  Carolina  continental  shelf.  Therefore,  it  may  be  a 
robust  feature,  not  specific  to  a  particular  model,  or  physics  pack¬ 
age,  or  hindcast.  Herein,  this  behavior  is  partially  attributed  to  the 
D1A  by  applying  the  same  hindcast  with  XNL.  In  this  study,  only 
one  model  (SWAN),  one  hindcast  (Lake  Michigan),  and  one  set  of 
physics  are  employed.  It  is  not  demonstrated  that  utilization  of 
XNL  in  another  model,  hindcast,  or  with  other  physics  would  yield 
the  same  improvement  to  model  agreement  with  observations. 
However,  it  is  noted  that  narrower  or  more  peaked  spectra  with 
XNL  is  a  robust  characteristic  of  comparisons  between  D1A  and 
XNL  (e,g,  Hasselmann  et  al.,  1985,  RW07,  Tolman,  2011),  so  it  is 
reasonable  to  expect  that  where  this  robust  bias  feature  exists, 
application  of  XNL  may  improve  agreement  with  observations. 

8.  Conclusions 

We  summarize  the  results  of  this  study  as  follows: 

•  From  analysis  of  hindcasts  for  Lake  Michigan  presented  in  this 
study,  as  well  as  prior  cited  works,  overprediction  of  energy 
below  the  spectral  peak— resulting  in  overly  broad  non-direc- 
tional  spectra— may  be  a  consistent  feature  of  3G  models  which 
utilize  the  Discrete  Interaction  Approximation  (D1A)  forS,,/4. 


•  Using  a  10-day  hindcast  for  Lake  Michigan  here,  and  in  prior 
idealized  computations  (Hasselmann  et  al.,  1985,  RW07),  we 
find  a  consistent  outcome:  the  models  which  use  exact  compu¬ 
tations  for  Sn/4  produce  more  narrow  frequency  spectra  than 
comparable  simulations  that  use  the  DIA  for  Sn/4. 

•  Four  methods  are  presented  for  quantifying  the  narrowness  of  a 
spectrum  in  frequency  space.  Those  denoted  as  Methods  C  and 
D  presented  here  are  found  to  be  most  useful,  insofar  as  they  are 
found  to  be  suitably  sensitive  to  the  narrowness  that  is  clear 
from  visual  inspection.  Method  C  is  a  frequency  analog  of  a 
method  proposed  by  Babanin  and  Soloviev  (1987,  1998a)  for 
quantifying  directional  width.  Method  D  employs  the  peaked¬ 
ness  factor  used  by  Goda  (1970,  1985). 

•  The  narrowness  quantity  associated  with  Method  C  agrees 
well— in  the  mean— with  observed  frequency  spectra  if  exact 
computations  are  used  forSn,4  in  the  10-day  hindcast  presented. 

•  Scatter  of  the  same  narrowness  quantity  for  the  same  simula¬ 
tion  is,  however,  worse  when  exact  computations  are  used  for 
Snl4* 
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Appendix  A.  Additional  error  statistics,  including  results  with 
alternate  physics 

Though  not  a  specific  goal  of  this  study,  it  is  useful  to  document 
the  performance  of  the  models  in  terms  of  conventional  quantities, 
wave  height  and  mean  period.  As  mentioned  earlier,  validation  of 
higher  order  moments  is  more  meaningful  in  cases  where  lower 
order  moments  are  in  good  agreement.  In  this  appendix,  we  in¬ 
clude  the  following  set  of  quantities: 

1.  Significant  waveheight,  Hm0  =  4^^  =  4y/J‘E(f)df. 

2.  Mean  spectral  period,  Tm0i,  =  fE(f)df/  JfE(f)df. 

3.  mean  spectral  period,  Tm,0t2  =  y/JE(f)df/  JPE(f)df. 

4.  Mean  spectral  period,  Tm_10  =  //  1  E(f)df/JE(f)df. 

5.  Spectral  narrowness  (or  peakedness)  parameters,  Q*.  Qb,  Qc.  and 
Qd  as  defined  in  Section  4. 

In  all  integrals,  the  lower  and  upper  bounds  on  integration  fmin 
and/ma*  are  implied. 

Statistics  used  to  quantify  error  are  the  following: 

1.  Bias,  i.e.  the  mean  error, 

2.  Normalized  bias,  the  bias  divided  by  the  mean  of  the 
observations. 

3.  Root-mean-square  error  (RMSE). 

4.  Slope  of  a  least-squares  fit  that  passes  through  the  origin. 

5.  Pearson’s  correlation  coefficient  (CC)  as  computed  by,  for  exam¬ 
ple,  Cardone  et  al,  (1996),  Ardhuin  et  al.  (2010), 
CC  -  ■ ,  <(Q;  W  M]i  -,  where  and  (  )  indicate  a  mean,  O  are 

Vi(0  0)2)V((M  M)2) 

observations  and  M  are  model  values. 

6.  Normalized  root-mean-square  error,  as  given  by  Ardhuin  et  al. 

(2010),  NRMSE  = 

7.  Scatter  index  SI,  the  standard  deviation  of  errors  divided  by  the 
mean  of  observations. 

It  is  also  useful  to  present  results  using  an  alternate  physics 
package,  especially  since  the  physics  used  in  this  paper  are  rela¬ 
tively  new.  This  is  done  here  using  the  Komen  et  al.  (1984)  physics 
(denoted  here  as  KHH)  with  the  minor  adjustment  suggested  by 
Rogers  et  al,  (2003).  [The  dependence  on  relative  wavenumber 
nk  *  2  is  used,  where  nk  =  1  is  favored  by  Komen  et  al.  (1984).  Here, 
Sdsif)  oc  (k/km)n\  and  km  is  the  mean  wavenumber].  The  physics 
have  deficiencies,  especially  with  regard  to  non-physical  interac¬ 
tion  between  windsea  and  swell,  and  non-physical  spectral  slope 
in  the  high-frequency  tail,  but  as  shown  by  Rogers  et  al.  (2003) 
and  RW07,  using  nk  -  2  can  be  highly  skillful  in  the  Lake  Michigan 
simulations.  For  these  physics,  calculations  were  performed  only 
with  the  DIA. 

For  sake  of  completeness,  we  present  statistics  with  and  with¬ 
out  the  sub-selection  of  spectra  according  to  eH  and  er  as  described 
in  Section  4. 

Statistics  for  the  three  types  of  model  runs  are  given  in  Table  2. 
The  results  for  the  DIA  and  XNl  based  runs  support  the  conclusions 
about  the  applicability  of  the  various  narrowness  parameters  pre¬ 
sented  in  Section  4.  The  results  for  the  DIA  based  runs,  using  the 
’old’  and  ‘new’  physics  package  are  close  to  each  other;  this  sug¬ 
gests  that  the  applicability  of  the  narrowness  parameters  does 
not  depend  on  the  type  of  physics  package.  The  results  also  show 
that  the  DIA-based  model  runs  have  better  performance  (in  terms 


of  integral  wave  height  and  period  measures)  than  the  XNL  based 
runs.  This  is  probably  due  to  the  fact  that  the  XNL  based  model  was 
not  calibrated,  but  we  deem  this  not  essential  for  our  purpose.  We 
also  present  statistics  for  a  longer  simulation,  that  used  by  RW07. 
These  are  given  in  Table  3,  For  this  74-day  simulation,  calculations 
were  performed  only  with  the  DIA.  As  already  concluded  for  the  re¬ 
sults  presented  in  Table  2,  the  type  of  input/dissipation  physics 
package  has  practically  no  effect  on  the  applicability  of  the  narrow¬ 
ness  parameters. 
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